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ABSTRACT 

We present a technique for bridging the gap between ideahzed inverse covariance weighted quadratic 
estimation of 21 cm power spectra and the real-world challenges presented universally by interferomet- 
ric observation. By carefully evaluating various estimators and adapting our techniques for large but 
incomplete data sets, we develop an optimal power spectrum estimation framework that preserves the 
so-called "EoR window" and keeps track of estimator errors and covariances. We apply our method 
to observations from the 32-tile prototype of the Murchinson Widefield Array to demonstrate the im- 
portance of a judicious analysis technique. Lastly, we apply our method to investigate the dependence 
of the clean EoR window on frequency — especially the frequency dependence of the so-called "wedge" 
feature — and establish upper limits on the power spectrum from z = 6.2 to z = 11.7. Our lowest limit 
is A(fc) < 0.3 Kelvin at 95% confidence at a comoving scale k = 0.046 Mpc~* and z = 9.5. 



1. INTRODUCTION 

In recent years, 21cm tomography has emerged as a 
promising probe of the Epoch of Reionization (EoR). As 
a direct measurement of the three-dimensional distribu- 
tion of neutral hydrogen at high redshift, the technique 
will allow detailed study of the complex astrophysical in- 
terplay between the intergalactic medium and the first 
luminous structures of our Universe. This will eventu- 
ally pave the way towards the use of 21 cm tomography to 
constrain cosmological parameters to exquisite precision, 
thanks to the enormi ty of the physical spac e within its 



reach (ple ase s e e, e.g., Furlanetto et al.'('2006');'Morales fc 



Wyithe( 2010| 

eTEol ( [201 3]) tor recent reviews) 



Pritchard & Locb^(2012};,Loeb &: Furlan- 
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To date, observational efforts have focused on measure- 
ments of the 21cm power spectrum. Such a measure- 
ment is exceedingly difficult. Sensitivity requirements 
are extreme, requiring thousands of hours of integ ration 
and large collecting areas (Mora les 2005; Bow man et al 
gOOg JLidz eTar i[2008l [Marker et al. 2010; Parsons et al 



|2012a r Adding to this challenge is the fact that raw 
sensitivity is insufficient — what counts is sensitivity to 
the cosmological signal above expected contaminants like 
galactic synchrotron radiation, which are three to four 
orders of magnitude brighter at the releva nt frequ encies 
jde Oliv eira-C osta et al.. 20081 pelic et al.|2008^ Bernardi 
efarp 009 Po beret al.|[Mfp . 

To deal with these challenges, numerous techniques 
have been proposed and implemented for foreground mit- 
igation and power spectrum estimatio n. These include 
foreground removal via parametric fits ( Wang et al.|2006[ 



Bowman et al.'2009 l|Liu et a l. 2009a'bp^ 

methods plarker et arpDUg| ',Chapmaii et al.l2012, 2013|), 
principal component analyses (Paciga et al. 2011; Liu &; 



principal component analyses (Paciga 
lTegmark'2012VMasui et al.'2013; Paciga ct al. 2015^111 
tering (Cleser ct al. 2008; Pctrovic & Oh 2011 [Parsons 
[et al.[[2012b[ ), treque ncy stac king (Cho ct al. 2012'), and 



quadratic methods (Liu fc Tegmark^^2011^ ^Dillon et al 



2013|[Shaw et al.|2013p . in almost all of these proposals 



foregrounds are separated from the cosmological signal 
by taking advantage of the differences in their spectra. 
Foregrounds are dominated by continuum processes and 
thus have smooth spectra. On the other hand, because 
the cosmological line-of-sight distance maps to the ob- 
served frequency of the redshifted 21 cm line, the rapid 
fluctuations in the brightness temperature distribution 
that are expected from theory will map to a measured 
cosmological signal with jagged, rapidly fluctuating spec- 
tra. When these spectral differences are considered in 
conjunction with instrumental characteristics, one can 
identify an "EoR window" : a region in Fourier space 
where power spectrum measure ments are expected to be 
relatively free from foregrounds (Datta et al.[[2010 Par- 
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Figure 1. The "EoR window," a region of Fourier space with rel- 
atively low noise and foregrounds, is thought to present the best 
opportunity for measuring the cosmological 21 cm power spectrum 
during the Epoch of Reionization. Here we show an example power 
spectrum from early MWA data, as a function of k±^ (Fourier 
components perpendicular to the line of sight) and fen (Fourier 
components parallel to the line of sight). More details on how 
we have calculated and plotted P{k^,k\\) are found in Section p| 
We schematically highlight the instrumental and foreground effects 
that that delimit the EoR window — the coldest part of this power 
spectrum. At low and high fcx, measurements are limited by an in- 
strument's ability to probe the largest and smallest angular scales, 
respectively. Limited spectral resolution causes similar effects at 
the highest fcii . As spectrally smooth sources, foregrounds inhabit 
primarily the low fcii regions. Thanks to chromatic instrumental 
effects, however, there is a slight encroachment of foregrounds to- 
wards higher fcii at higher k^, in what has been colloquially termed 
the "wedge" feature. 



sons et aL " 2012b: Vcdantham ct al."2012: Morales et aL| 
2012| [Trott ct al. ^ ,2012, ,Thyagarajan ct al.„2013J . This 
is shown schematically in l^'igurejIJ where we have used 
early Murchison Widefield Array (MWA) data to esti- 
mate the power spectrum as a function of k±^ (Fourier 
mode perpendicular to the line-of-sight) and fcu (Fourier 
mode parallel to the line-of-sight). More details regard- 
ing this figuew are provided in Section |3] for now we 
simply wish to draw attention to the existence of a rela- 
tively contaminant-free region in the middle of the fc_L-fcj| 
plane. This clean region is what we denote the EoR win- 
dow. 

The EoR window is generally considered the sweet spot 
for an initial detection of the cosmological 21 cm power 
spectrum, and constraints are likely to degrade away 
from the window. At high k^_ (i.e., the finest angular 
features on the sky), errors increase due to the angular 



resolution limitations of one's instrument. For an inter- 
ferometer, this resolution is roughly set by the length of 
the longest baseline. Conversely, the shortest baselines 
define the largest modes that are observable by the in- 
strument. Errors therefore also increase at the lowest fcj^ 
where again there are few baselines. 

A similar limitation defines the boundary of the EoR 
window at high /cy. Since the spectral nature of 21cm 
measurements mean that different observed frequencies 
map to different redshifts, the highest fcy modes are in- 
accessible due to the limited spectral resolution of one's 
instrument. At low k\\, one probes spectrally smooth 
modes — precisely those that are expected to be fore- 
ground contaminated. Thus there is another boundary 
to the EoR window at low fcy . 

A final delineation of the EoR window is provided by 
the region labeled as the "wedge" in Figure[T] The wedge 
feature is a result of an interplay between angular and 
spectral effects. Simulations have shown that the wedge 
is the effect of chromaticity in one's synthesized beam 
(which is inevitable when an interferometer is used to 
survey the sky). This chromaticity imprints unsmooth 
spectral features on measured foregrounds, resulting in 
foreground contamination beyond the lowest k\\ modes 
even if the foregrounds themselves are spectrally smooth. 
Luckily, this sort of additional contamination follows a 
reasonably predictable pattern in the fc_L-fc|| plane, and in 
the limit of intrinsically smooth foregrounds, the wedge 
can be shown to extend no farther than the line 



fell 



sn\9 



field- 



Dm{z)E{z) 
Dh{1 + z) 



(1) 



where Dh = c/Hq, E{z) = ^n^{l + zf + nf^, 

Dm{z) = /q dz'/E{z'), ^fioid is angular radius of the the 
field-of-vie w, and c. Ho , flrr,., and ft\ have their usual - 



meanings (iDatta et al. 2010 [Vedantham et al 



Morales etai.||2012[ [Trott et al.||2012p . Intuitive 



2012 

_^_^ tHe 

toreground-contaminated wedge extends to higher fcy at 
higher k± because the high k± modes are probed by the 
longer baselines of an interferometer array, which have 
higher fringe rates that more effectively imprint spectral 
structure in the measured signals. For an alternate but 
equivalent explanation in terms of delay mode s , pleas e 
see the illuminating discussion in Parsons et al. (2012b). 
The concept of an EoR window is important m that it 
provides relatively strict boundaries that separate fairly 
foreground-free regions of Fourier space from heavily 
foreground-contaminated ones. It therefore provides one 
with the option of practicing foreground avoidance rather 
than foreground subtraction. If it turns out that fore- 
grounds cannot be modeled well enough to be directly 
subtracted with the level of precision required to detect 
the cosmological signal, foreground avoidance becomes 
an important alternative, in that the only way to robustly 
suppress foregrounds is to preferentially make measure- 
ments within the EoR window. Likely, some combina- 
tion of the two strategies — foreground subtraction and 
foreground avoidance — will prove useful for the detec- 
tion of the 21 cm power spectrum. Of course, measure- 
ments within the EoR window are still contaminated by 
instrumental noise, but fortunately the noise integrates 
down with further observation time (as long as calibra- 
tion errors and other instrumental systematics can be 
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sufficiently minimized). Observationally, it is encourag- 
ing that the EoR window has now been shown to be free 
of fore grounds to better than one part in a hundred in 
power ( (Pober et alT]|2013| . 

As experimental sensitivities increase, however, one 
must take care to preserve the cleanliness of the EoR 
window to an even higher dynamic range. There are 
several ways in which our notion of the EoR window 
may be compromised. First, as experiments integrate 
in time and acquire greater sensitivity, we may discover 
that our approximation of spectrally smooth foregrounds 
is insufficiently good for a detection of the (faint) cosmo- 
logical signal. In other words, foreground sources may 
have small but non-negligible high fcy components in their 
spectra that have thus far gone undetected. This would 
translate into a smaller-than-expected EoR window. In 
addition, even intrinsically smooth foregrounds may ap- 
pear jagged in a real measurement because of instru- 
mental effects such as imperfect calibration. The precise 
interferometer layout may also result in unsmooth arti- 
facts that arise from combining d ata from non-redundant 
baselines (Hazelton et al. 2013). Finally, suppose that 
the aforemientioned ettects are negligible and that the as- 
sumption of spectrally smooth foreground emission con- 
tinues to hold. The EoR window still cannot be taken 
for granted because non-optimal data analysis techniques 
may result in unwanted foreground artifacts in the re- 
gion. For the EoR window to exist at all, it is essential 
that power spectra are estimated in a rigorous fashion, 
with well-understood statistics. 

The goal of this paper is to minimize unwanted data 
analysis artifacts by establishing methods for power spec- 
trum estimation that are both robust and optimal. Pre- 
vious efforts have rarely met both criteria: either the 
methods are robustly applicable to data with real-world 
artifacts but fail to achieve optimal (or even rigorously 
computable) error properties, or provide an optimal 
framework but ignore real-world complications. In this 
paper we ex tend the rigorous framework described in |Liu| 
fc Tegmark] (|2011[ ) and [Dillon et a l.| (|2013[) to deal with 
reai-world etiects. The result is a computationally fea- 
sible approach to analyzing real data that not only pre- 
serves the cleanliness of the EoR window, but also rigor- 
ously keeps track of all relevant error statistics. 

To demonstrate the applicability of our approach, we 
apply our techniques to early data from the Murchison 
Widefield Array (MWA). These data were derived from 
~ 22 hours of tracked observations using an early, 32- 
element prototype array. The results are therefore not 
designed to be cosmologically competitive, but instead 
illustrate the rigor that will be required for an eventual 
detection of the EoR while also providing new measure- 
ments on the "wedge" feature that delineates the EoR 
window. 

This paper is organized as follows. In Section[2]we dis- 
cuss various real-world obstacles that must be dealt with 
when analyzing real data, and how one can overcome 
them while maintaining statistical rigor. We then apply 
our methods to MWA data in Section [3] as a "worked ex- 
ample" , highlighting the importance of various subtleties 
of power spectrum estimation. In Section HI we present 
some results from the data, emphasizing the agreement 
between theoretical expectations and our observations 
of the foreground wedge (particularly regarding the fre- 



quency dependence of the wedge) . We also present upper 
limits on the cosmological 21 cm power spectrum over the 
broad redshift range of2: = 6.2toz = ll.l. Finally, we 
summarize our conclusions in Section [H 

2. SYSTEMATIC METHODS FOR DEALING WITH 
REAL- WORLD OBSTACLES 

To understand the gap between an analysis framework 
for idealized observations and any real- world data set, we 
enumerate and address six different obstacles that rather 
universally affect real data. Our goal in this section is to 
meet the challenges presented by these obstacles while 
maintaining the advantages of the optimal framework, 
which we reiterate in Section |2.1[ especially the abil- 
ity to minimize and precisely quantify the uncertainties 
in the measurements. In the following sections, we ad- 
dres s the problems presented by large data volumes (Sec- 
tion 2.2), uncertainties in th e pr operties of contaminants 
such as f oreg rounds (Section 2.3), incomplete uv coverage 
(Section |2.4[) , radio frequency interference (RFI) flagging 
(Section J2J)j ) , foreground leakage into the EoR window 
(Section |2.6[ ), an d bi nning to spherically averaged power 
spectra (Section 2.7). 



2.1. A systematic framework for analyzing idealized 
observations 



In this section , we briefly review the formalism of |Liu fc] 
Tegmark (2011) for optimal power spectrum estimation, 
which was adapted for 21 cm tomography from similar 
techniques used in galaxy survey and cosmic microwave 
background ana lysis (Tegmark 1997b| [Bond et alT 1998 
Tegmark et al.[[l998[ [2004p. lor now, wc do not include 



real-world ettects such as missing data from RFI flag- 
ging, and the purpose of later sections is to extend the 
formalism to take into account these complications. 

In 21cm tomography, one typically wishes to measure 
both the spherically-binned power spectrum Psph(^), de- 
fined by 

(f*(k)f (k')) ^ (27r)3p,ph(fc)5(k - k'), (2) 

and the cylindrically-binned power spectrum 
Pcyi{k±,k\\), defined by 

(f*(k)r(k')) = (27r)3p,yi(fci,fc||)5(k-k'), (3) 

with r(k) signifying the spatial Fourier transform of the 
21 cm brightness temperature field T(r), k denoting the 
spatial wavevector with magnitude fc, and components 
k± and ku as the components perpendicular and paral- 
lel to the line-of-sight, respectively. The angled brackets 
(• • • ) represent an ensemble average. The spherical power 
spectrum is useful for comparing to theoretical models, 
since it is obtained by angularly averaging over spher- 
ical shells in Fourier space, and thus makes the cosmo- 
logically relevant assumption of isotropy. The cylindrical 
power spectrum is useful for identifying instrumental and 
foreground effects, which possess a cylindrical symmetry 
rather than a spherical one. Typically, the cylindrical 
power spectrum is produced first as a tool for foreground 
isolation (i.e., to identify the EoR window), and then 
subsequently binned into a spherical power spectrum. 
This section concerns the estimation of the cylindrical 
power spectrum. Optimal binning techniques to go from 
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the cylindrical spec trum to the spherical spectrum are 
discussed in Section [2?7l 

In estimating a power spectrum from data, one must 
necessarily discretize the problem. We make the approx- 
imation that the power spectra are piecewise constant 
functions, such that we can describe them in terms of a 
vector of bandpowers with components p"' , where 



Pa = -^ cyl(rS_Li "^11 )■ 



(4) 



It is the bandpowers and their error properties that one 
wishes to estimate from the data, which come in the form 
of a data vector x. Intuitively, one can think of the data 
vector as a list of the 21 cm brightness temperatures mea- 
sured at various locations in a three-dimensional "data 
cube". Rigorously, we define each element of the data 
vector (i.e., each voxel of the data cube) as 



X, = y"r(r)V^,(r)d3r, 



(5) 



with ■(/'i(r) being the pixelization kernel and T(r) as the 
(continuous) three-dimensional 21 cm brightness temper- 
ature fiel(p] In this paper we take the i*'' pixelization 
kernel V'i(iv to be a boxcar function centered on the i*'* 
voxel of the data. 

To estimate the a"* bandpower from the data vector, 
we first form a quadratic estimator of the form 



g„ =i(x - m)*C-iC,„C-i(x - m) 



(6) 



where m = (x) is the mean of the data, C = (xx*) — 
(x)(x)* is its covariance, Cjun^ is the component of the 
covariance "junk" /contaminants (to be defined in the fol- 
lowing section), and C_q is the derivative of the covari- 
ance with respect to the a*'' bandpower. Since we are ap- 
proximating the power spectrum as piecewise constant, 
we have 



-^junk 



E 



PaC,, 



(7) 



Combined with Equation (l5| , this expression can be used 
to derive explicit forms for C_q, which reveals that the 
matrix essentially Fourier transforms and bins the data 



( Liu fc Tegmark||2011 
C^Q can be thought o 



2013). Intuitively, 
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as the response in the data co- 
variance C to the bandpower pa. Thus, as long as one 
selects an appropriate form for C_q,, the formalism of this 
section can also be used to directly measure the spherical 
power spectrum. However, as we discussed above, in this 
paper we choose to first estimate the cylindrical power 
spectrum as an intermediate diagnostic step, to quantify 
and mitigate foregrounds better. 

^^ Of course, instrumental noise and foregrounds do not prop- 
erly reside in a cosmological three-dimensional volume: noise is 
introduced in the electronics of the system, whereas foregrounds 
are "nearby" and only appear in the same location in the data 
cube as our cosmological signal by virtue of their frequency de- 
pendence. However, there is a gain in convenience and no loss of 
generality in assigning a noise and foreground contribution to each 
voxel, pretending that those contaminants also live in the observed 
cosmological volume. 



Once the g^s have been formed, they need to be nor- 
malized using a suitable invertible matrix M to form the 
final bandpower estimates: 



P = Mq, 



(8) 



where we have grouped the bandpower estimates Pa into 
a vector p (and similarly grouped the coefficients qa and 
q) , with the hat ( ^ ) signifying the fact that we have 
formed an estimator of the true bandpowers^ We shall 
discuss different choices of M in Section 12.61 

To understand the uncertainty in our estimates, we 
compute several error properties. The first is the covari- 
ance matrix of the final measured bandpowers: 



S=(pp*)-(p)(p)* = MFM*, 



(9) 



where we have introduced the Fisher matrix F, which 
has components 



F^p = itr[C-iC,„C-iC,^]. 



(10) 



The Fisher matrix also allows us to relate our estimated 
bandpowers p to the true bandpowers p via the window 
function matrix W: 



(p) = Wp, 
where W can be shown to take the form 
W = MF. 



(11) 



(12) 



If we choose M such that the rows of W each sum to 



unity. Equation (11) shows that each bandpower esti- 
mate can be thought of as a weighted average of the 
truth, with weights given by each row (each window func- 
tion). Even with this normalization requirement, there 
are still many choices for M. We discuss the various 
options and tradeoffs in Section [2^1 

Whatever the choice of M, our estimator has optimal 
error properties in the sense that if p in Equation ( 11 ) is 



used to constrain parameters in some theoretical model, 
those measured parameters will have t he smallest p ossi- 
ble error bars given the observed data ( Tegmark"|T997bJ . 
Our goal in the following sections will be to ensure that 
both these small error bars and our ability to rigorously 
compute them are preserved in the face of real-world dif- 
ficulties. 

2.2. A real-world obstacle: data volume 

Perhaps the most glaring difficulty presented by the 
ideal technique outlined above is its computational cost. 
Much of that cost arises from the inversion of the data co- 
variance matrix C in Equations (rol and ( |10p , in addition 
to the multiplication of C and matrices of the same size. 
Both of these operations scale like 0{N^), where N is the 
number of voxels in each data vector. The computational 
cost makes taking full advantage of current generational 
interferometric data prohibitive, not to mention upcom- 
ing observational efforts that expect to produce 10^ or 
more voxels of data. 

^® Note that q, p, and M live in a different vector space than 
X, C, and C,c«. The former are in a vector space where each 
component refers to a different bandpower, whereas the latter are 
in one where different components refer to different voxels. 
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One would like to retain the information theoretic ad- 
vantages of the quadratic estimator method and its abil- 
ity to precisely model errors and window functions, with- 
out 0{N^) complexity. The s olution to this prob lem, de- 
veloped and demonstrated in Dillon et al. (2013), comes 
from taking advantage of a number ot symmetries and 
approximate symmetries of the survey geometry and the 
covariance matrix, C, and can accelerate the technique 
to O(A^logiV). 

The fast method relies on assembling the data into a 
data cube with rectilinear voxels amenable to manipula- 
tion with the Fast Fourier Transform. This is equivalent 
to the assertion that each voxel represents an equal vol- 
ume of comoving space, an approximation that relies on 
two restrictions on the data cube geometry. First, the 
range of frequencies considered must be small enough 
that Dc{z) (the line-of-sight comoving distance, equal to 
Dm{z) above in a spatially flat universe) is linear with 
V. Generally, one should limit oneself to analyzing the 
power spectrum of redshift ranges short enough that the 
evolution of the power spectrum durin g reionization can 
be neglected. This range, suggested by Mao et al. (20081 
to be Az < 0.5, makes the approximation ot a linear re- 
lationship between v and D(.{z) better than one part in 
10'^ at the redshifts of interest to 21 cm cosmology. 

Second, the assumption of equal volume voxels relies 
on the flat sky approximation. To achieve this the area 
surveyed can be broken into a number of subfields, each 
a few degrees on a side, for which the curvature of the 
sky can be neglected. As long as the angular extent of 
the data cube is smaller than ^ 10°, the flat sky approx- 
imation is correct to a few parts in 10"^. 

By analyzing a rectilinear volume of the universe, all 
steps in calculating the band powers qa can be performed 
quickly by exploiting various symmetries and taking ad- 
vantage of the Fast Fourier Transform. The model for C 
can be broken up into a number of independent matrices 
representing signal, n oise, and foregrounds. E ach of these 
models, developed by Liu & Tegmark (2011), is well ap 



proximated by a sparse matrix m a convenient com bina- 
tion of real and Fourier spaces ( Dillon et al.||2013 ). As a 
result, multiplica tion of a vector by C can be performed 
in 0{N\ogN). [Dillon et al] \2Q12,\ showed how that 
speed-up can be parlayed into a method for quickly cal- 
culating Qa using the Conjugate Gradient Method. The 
rapid convergence of the iterative method for calculat- 
ing C^^x can be ensured by the application of a pre- 
conditioner which relies on the spectral smoothness of 
foregrounds and the fac t that they are well d escribed by 
only a few eigenmodes ( Liu k, Tegmark|2012 |. Then, by 
randomly simulating many data vectors trom the covari- 
ance C and calculating qa from each, the Fisher matrix 
can be estimated from the fact that 

F-(qq*)-(q)(q*), (13) 

which follows from Equation (|9|. All of this together 
allows for fast, optimal power spectrum estimation — 
including error bars and window functions — despite the 
challenge presented by an enormous volume of data. 

2.3. A real-world obstacle: uncertain contaminant 
properties 

If one had perfect knowledge of the foreground contam- 
ination in the data cube, the problem of foreground con- 



tamination would be trivial; one would simply perform a 
direct subtraction of the foregrounds from the data vec- 
tor X. Unfortunately, our knowledge of foregrounds is 
far from perfect, particularly at the level of precision re- 
quired for a direct detection of the cosmological 21 cm 
signal. Because of this, the estimator shown in Equation 
([6]) in fact combines several different foreground subtrac- 
tion steps in an attempt to achieve the lowest possible 
level of foreground contamination: 

1. A direct subtraction of a foreground model from 
the data vector. This is given by x—m. To see this, 
note that the data vector can be thought of as being 
comprised of the cosmological 21 cm signal X21, the 
foregrounds Xfg, and the instrumental noise n. On 
the other hand, the mean data vector 



m = W = (X21) 



(Xfg) + (n) = (xfg). (14) 

contains only the foreground contribution, because 
we are interested in the fluctuations of the 21 cm 
signal, so the cosmological signal has zero mean, 
as does the instrumental noise (in the absence of 
instrumental systematics). Note that because the 
mean here is the mean in the ensemble average 
sense (as opposed to just the spatial mean), m rep- 
resents a full spatial and spectral model of the fore- 
grounds. 

2. Since the foregrounds also appear in the covari- 
ance matrix, the action of C~^ is to downweight 
foreground-contaminated modes, exploiting fore- 
ground properties such as smooth frequency depen- 
dence. 



3. Subtracting the term itr[CjunkC ^C^qC 
inates the bias from contaminants. 



elini- 



4. Finally, the binning of the cylindrical power spec- 
trum to the spherical power spectrum provides 
yet more foreground suppression. Foregrounds are 
distributed in select regions on the fcj_-A:|| plane 
(i.e., outside the EoR window) in patterns that do 

not lie along contours of constant k = 1 /fc^ + k?,. 

Thus, when binning along such contours to produce 
a spherical power spectrum, one can selectively 
downweight parts of the contour with greater fore- 
ground contamination, which constitutes a form of 
foreground cleaning. Roughly speaking, this corre- 
sponds to taking advantage of the fact that fore- 
grounds have a cylindrical symmetry in Fourier 
space, w hereas the signa l is spherically isotropic 
(Morales fc Hewitt|[2004| . We do note, however, 
that the formalism we introduce in Section 12.71 is 
general enough to use any geometric differences be- 
tween foregrounds and signal. 

Of these foreground mitigation strategies, the first and 
third are direct subtractions (in amplitude and power, 
respectively), whereas the second and the fourth act 
through weightings. The former group represent opera- 
tions that are particularly vulnerable to incorrectly mod- 
eled foregrounds. To see this, recall that the foregrounds 
are expected to be larger than the c osmologi cal signal 
by three or four orders of magnitude (de Oliveira-Costa] 



DILLON, LIU, WILLIAMS, ET AL. 



eFaL]|2008{ |Jelic et aL|[2008l |Bernardi et ar][2009l |Pober| 



et al. 12013 ). Thus, when performing direct subtractions, 
low-level, unaccounted-for inaccuracies in the foreground 
model can translate into extremely large biases in the fi- 
nal results. In addition, significant numerical errors may 
arise from the subtraction of two large numbers (the data 
and the foregrounds) to obtain a small number (the mea- 
sured cosmological signal). 

Our goal for the rest of the section is to immunize our- 
selves against biases from direct subtractions. Of the 
direct subtraction steps list above, the Step 1 is likely to 
be relatively harmless for two reasons. First, it is imme- 
diately followed by the C^^ downweighting. The down- 
weighting mitigates the effects of inaccuracies in model- 
ing, for the C~^ tends to gives less weight to precisely 
the modes that have the largest foreground amplitudes, 
and therefore would be the most susceptible to modeling 
errors in the first place. In addition, the uncertainty in 
foreground properties in those regions of the k±-k^\ plane 
result in large error bars there, providing a convenient 
marker of the untrustworthy parts of the plane, effec- 
tively demarcating the boundaries of the EoR window. 
For these two reasons. Step 1 is unlikely to be an issue, 
at least not inside the EoR window. 

More worrisome is Step 3, where the power spectrum 
bias of contaminants is subtracted off. If we define "con- 
taminants" to be "everything but the cosmological 21 cm 
signal", there are two potential sources of bias: fore- 
grounds and noise. The subtraction of these biases is 
not followed by a downweighting analogous to the appli- 
cation C~^ in Step 1. Moreover, whereas one could argue 
that the foreground bias is likely to be large only outside 
the EoR window, the noise bias will spread throughout 
the fcj_-fcj| plane. This noise bias will also be quite large, 
as current experiments are firmly in the regime where 
the signal-to-noise is below unity. It would therefore be 
advantageous to avoid bias subtractions altogether if pos- 
sible. 

To avoid having to subtract foreground bias, we sim- 
ply redefine what we mean by contaminants/junk. If 
we modify our mission to be one where we are measur- 
ing the power spectrum of total sky emission instead of 
the power spectrum of the cosmological 21 cm signal, the 
foreground contribution to the bias term no longer ex- 
ists, as foregrounds now count as part of the signal we 
wish to measure. Of course, nothing has really changed, 
for we have simply ignored the subtraction of the fore- 
ground bias by redefining what we mean by "contami- 
nants". Within the EoR window, this should result in 
little degradation of our final constraints, for in this re- 
gion foreground contamination is expected to be negligi- 
ble, and the power spectrum of the cosmological signal 
should be essentially identical to the power spectrum of 
total sky emission. In any case, this is an assumption 
that can be checked in the final results, and represents a 
conservative assumption throughout Fourier space since 
foreground power is necessarily positive. 

In contrast, escaping to the safe confines of the EoR 
window alone is not sufficient to eliminate the instrumen- 
tal noise portion of the bias term, for the instrumental 
noise bias pervades the entire k±-kH plane. To elimi- 
nate the noise bias, one can choose to compute not the 
auto-power spectrum of a single data cube with itself, 
but instead to compute the cross-power spectrum of two 



data cubes that are formed from data from interleaved 
(i.e., odd and even) time samples. Since the instrumen- 
tal noise is uncorrelated in time, this has the effect of 
automatically removing the instrumental noise bia^^ 

More explicitly, we can form a bandpower estimate of 
the cross-power spectrum by simply computing 



^ross 
Pa 



:X1E"X2, 



(15) 



where xi and X2 are the data vectors for the two time 
inter-leaved data cubes, and for notational brevity we 
have defined E" = i X^/j ^■^a/3C-iC,/3C^i. To see that 
the cross-power spectrum has no noise bias, let us decom- 
pose the data vectors Xj into the sum of s and n^, the 
signal and noise components respectively, where the sig- 
nal component has no index because it does not vary in 
time (note also that following the discussion above, any 
true sky emission counts as signal, so that s = X21 -|-Xfg). 
Inserting this decomposition into the preceding equation 
and taking the expectation value of the result gives 



(p^J°-)-((s + ni)*E"(s + n2)) 
= (s*E"s) + (ni)*E"s 
+ s*E"(n2) + (niE"n2) 



= (s*E"s), 



(16) 



where the last equality holds because the instrumen- 
tal noise has zero mean, i.e. (rii) = 0, and no cross- 
correlation between different times, i.e. (nin2) — 0. The 
resulting estimator depends only on the power spectrum 
of the signal, and there is no additive bias. 

Importantly, however, we emphasize that while we 
have eliminated noise bias by computing a cross-power 
spectrum, we have not eliminated noise variance. In 
other words, the instrumental nosie will still contribute 
to the error bars. To see this, consider the variance in 
our estimator, which is given by 



^cross /;;:x;ross;;rcross\ 

^aP —\Pa P/3 I 



(P'a"")(p'r") 



(x*E"x2X*E'5x2) - (x* E"X2) (x* E-^xs) (17) 



The second term simplifies to 



{p. 



;C:cross\ /x::^cross 



)=^(xlxi)(xJx^)E 

ijkl 



ij^kl- 



(18) 



^^ The reader may object to this by (correctly) pointing out that 
there exist errors that are correlated in time, with calibration er- 
rors being a prime example. The result would be a cross-power 
spectrum that still retained a bias. However, this does not inval- 
idate the cross-power spectrum approach, in the following sense. 
While biases will make our estimates of the power spectrum imper- 
fect, these estimate will not be incorrect — the final (biased) power 
spectra will still represent perfectly rigorous upper limits on the 
cosmological power, provided we are conservative about how we 
estimate our error bars. We will discuss how to make such conser- 
vative error estimates later on in this section and in Section [sTsl 
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Similarly, the first term is equal to 

ijkl 



= ^((xlxi)(x^x?,) + (xlx5=)(x^2x: 

ijkl 



+ (xlx^)(x{x^)JEf^.Ef,, (19) 

where in the last equality we assumed Gaussian dis- 
tributed data to simplify the four-point correlation. Our 
bandpower covariance is now 



ijkl 



(xl4)(xix^))Ej;.Ef,. 



(20) 



The first term in this expression consists only of auto- 
correlations, which contain both noise and signal: 

(xix*i) = ((s+ni)(s*+n*))-(s)(s)* = S + N = C, (21) 

where we have define d C to be the total data covariance 
(as defined in Section 2.1 ), S = (ss*) — (s)(s)* is the sky 
signal covariance (as per the discussion earlier in this 
section), and N = (riin^) = (n2n2) is the instrumental 
noise covariance. We have assumed that there is no cor- 
relatiorpj between the sky emission and the instrumental 
noise, so that (sn*) = (sria) = 0. 

The second term in our bandpower covariance consists 
only of cross-correlations, and thus contains no noise co- 
variance: 



(xiX*) = ((s + ni)(s*-|-n*))-S. 



(22) 



Putting everything together, we obtain 

S^J*^ = tr[CE"CE'3] + tr[SE"SE'^] . (23) 

This, then, is the error covariance of our cross power 
spectrum estimator. It gives less variance than the ex- 
pression for the auto power spectrum, which in the no- 
tation of this section takes the form 



S^"^'° = 2tr[CE"CE''] 



(24) 



Despite this difference between equations [23] and [24j one 
may conservatively opt to use the above covariance ma- 
trix for the auto-power spec trum to estimate error bars 
even when using Equation ( 15 ) to estimate the power 
spectrum itself. In fact, it may be prudent to make 
this choice because there exists the possibility that the 
noise between interleaved time samples may not be truly 
uncorrelated, making the true errors closer to those de- 
scribed by S**" °. In our worked example with MWA 
data in Section [31 we will conservatively use Equation 
(24) to estimate the errors of our cross-power spectrum. 
In summary, uncertainties in noise and foreground 
properties make it desirable to avoid trying to extract 

^^ Note that this assumption has nothing to do with whether or 
not the instrument is sky-noise dominated. A sky-noise dominated 
instrument will have instrumental noise whose amplitude depends 
on the sky temperature, but the actual noise fluctuations will still 
be uncorrelated with the sky signal. 



weak signals by performing subtractions between two 
large numbers (the contamination-dominated data and 
the possibly inaccurate contaminant models). Mathe- 
matically, the greatest concern comes with the subtrac- 
tion of the noise and foreground biases from power spec- 
tra estimates. To deal with the residual noise bias, one 
may evaluate cross-power spectra between interleaved 
time samples rather than auto-power spectra. To deal 
with the foreground bias, one can conservatively elect 
to simply leave it in when placing upper limits on the 
cosmological signal, and rely on the robustness of the 
EoR window to separate out the foregrounds from the 
cosmological 21 cm signal. In effect, one can practice 
foreground avoidance rather than foreground subtrac- 
tion, since the former (if it is sufficient for a detection 
of the cosmological signal) will be more robust than the 
latter in the face of foreground uncertainties. Finally, 
as a brute-force safeguard, to quantify such uncertain- 
ties, one can always vary the foreground mod el us ed in 
power spectrum estimation, as we do in Section 3.3 when 
we apply our methods to the worked example of MWA 
data. 

2.4. A real-world obstacle: incomplete uv-coverage 

While the methods of the previous section allow one to 
alleviate the effects of foreground modeling uncertainty, 
it is impossible to avoid the fact that real interferome- 
ters are imperfect imaging instruments. This is because a 
real interferometer will inevitably have uw-coverage that 
is non-ideal in two ways. First, the coverage is non- 
uniform, resulting in images that have been convolved 
with non-trivial synthesized beam kernels. Second, the 
uu-coverage is incomplete, in that certain parts of the 
uu-plane are not sampled at all. The idealized methods 
of Section |2.1| deals with neither problem, and in this 
section with augment the formalism to rectify this. 

Assume for a moment that uv coverage is complete (so 
that there are no "holes" in the uu-plane), but not nec- 
essarily uniform. In such a scenario, one has measured 
an unevenly weighted sample of the Fourier modes of 
the sky. The effect of this non-trivial weighting needs to 
be accounted for when measuring the power spectrum, 
since uv coordinates roughly map to fcj^. A failure to 
do so would therefore result in the final power spectrum 
estimate being multiplied by some function of fcj^ corre- 
sponding to the uv distribution. 

Put another way, the uv distribution of an interferom- 
eter defines its synthesized beam, the kernel with which 
the true sky has been convolved in the produ ctio n of our 
image data cube. The equations of Section |2.1| assume 
that this convolution has already been undone. Thus, we 
must first perform this step, which in our notation may 
be written as 

x^B-^x', (25) 

where x' represents the convolved data vector, B is the 
convolution matrix encoding the effects of the synthe- 
sized beam, and x is the processed data vector that is 
fed into Equation ([6]). Note that this application of B~^ 
is meant to undo only the effects of the synthesized beam, 
not the primary beam. 

The above method assumes that the matrix B is in- 
vertible. In practice, this will likely not be the case as 
parts of the uv plane will be missed by the interferon!- 
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eter, resulting in a singular B matrix. In what follows, 
we will present two different ways to deal wi th th is. The 
first is to modify the equations of Section |2.1| so that 
they accept the convolved images (the "dirty maps") as 
input. Since all the statistical information relevant to the 
power spectrum are encoded in the covariance matrix, we 
simply have to make the replacement 



C^(xx*)-(x)(x>* 
This amounts to 



B((> 



(x')(x'; 



B' =BCB* 



(26) 



(27) 



Of course, changing the covariance matrix also changes 
C_Q, and we must propagate this change. Differentiating 
the preceding equation with respect to the bandpower p^ 
gives the substitution 



^.a 



BC ckB 



(28) 



Since C^q is the response of the data covariance C to the 
bandpower pa , this is simply a statement of the fact that 
if our data consists of dirty maps, the revised C.q. matrix 
should encode the response of a dirty map's data covari- 
ance to the b and power . W ith the substitutions given 



by Equ ation s ( 27 1 and ( 28 1 , the rest of the equations of 
Section 2.1| can be used unchanged. In the limit of an 
invertible B matrix, it is straightforward to show that 
this is equivalent to using Equation (25). 

The second method for dealing with a singular B, 
which was proposed in Dillon et al. (2013), is to replace 
the ill-defined inverse matrix B~^ with a pseudoinverse 
given by 

(29) 



n(B + 7UU^) n. 



where 7 is a non-zero but otherwise arbitrary real num- 
ber, and n is a projection matrix given by 



n = i-u(u^u)-iu 



(30) 



The matrix U specifies which modes on the sky are miss- 
ing in the data as a result of unobserved pixels on the 
uw-plane. It is constructed by computing the responses 
(on the sky) of each unobserved uv pixel individually and 
storing each response as a column of U. As an example, 
in the fiat-sky approximation the U matrix would have 
a sinusoid in each column, corresponding to the fringes 
that would have been observed by the interferometer had 
data not been missing in a particular uv pixel. If these 
modes were present in the covariance model (which might 
be the case, for example, if the covariance were con- 
structed by modeling data from a different interferometer 
with different uv coverage), then the inverse covariance 
C""'^ in our estimator needs to be similarly replaced with 
the pseudoinverse: 



n (c + 7UU 



n. 



(31) 



Importantly, the pseudoinverse can be quickly multiplied 
by a vector using the previously discussed conjugate gra- 
dient method. Its usage therefore does not sacri fice any 
of the speedups that were identified in Section |2.2| for 
dealing with large data volumes. 

2.5. A real-world obstacle: missing data from RFI 



In any practical observation, the presence of narrow- 
band RFI will mean that certain RFI-contaminated fre- 
quency channels will need to be flagged as outliers and 
omitted from a final power spectrum analysis. The re- 
sult, once again, is the presence of gaps in the data, only 
this time the missing modes are complete frequency chan- 
nels. However, the pseudoinverse formalism of the pre- 
vious section is quite flexible in that modes of any form 
can be projected out of the analysis. Thus, to correctly 
account for RFI-flagged data, one simply uses the pseu- 
doinverse in exactly the same way as one does to account 
for missing uv data. 

2.6. A real-world obstacle: foreground leakage into the 
EoR window 



As Equation ( 11 ) showed, estimates of the power spec- 
trum are not truly local, in the sense that every band- 
power estimate Pa corresponds to a weighted average 
of the true power spec trum, with weigh t s spe cifled by 



the window functions. Liu & Tegmark (2011) showed 



that these window functions can be quite broad, par 
ticularly in regions with high foreground contamination. 
There is thus the danger that foreground power could 
leak into the EoR window. Because the foregrounds are 
so much brighter than the cosmological signal, even a 
small amount of leakage could compromise the cleanli- 
ness of the EoR window. 

Fortunately, one can exert some control over the shape 
of the window functionfF^by making wise choices regard- 
ing the form of M in Equation ([8]), which in turn gives 
the window functions via W = MF. As discussed above, 
M must be chosen such that the rows of W sum to unity. 
Beyond that requirement, however, an infinite number of 
choices are admissible. One choice would be M = F~^, 
which gives W = I (i.e., delta function windows). This 
would certainly minimize the amount of leakage into the 
EoR window, but it comes at a high price: the result- 
ing error bars on the power spectrum measurement — the 
diagonal elements of S from Equation ^ — tend to be 
large, refiecting the data's inability to make highly lo- 
calized measurements in Fourier space when the survey 
volume is finite. 

On the other extreme, the error bars predicted by S 
can be shown to be their smallest poss ible if M is taken 
to be diagonal ( Tegmark et al.| 2002). However, this 
gives broader window functions, for it is via the smooth- 
ing/binning effect of these broad window functions that 
the small errors can be achieved. One can also argue 
that the level of smoothing dictated by this approach 
is excessive, since the resulting bandpowers have posi- 
tively correlated errors. (To see this, note that up to a 
row-dependent normalization, the error covariance ma- 
trix takes the form S ^ F. Since all elements of a 
Fisher matrix must necessarily be non-negative, this im- 
plies that all cross-covariances of the estimated bandpow- 
ers have positively correlated errors unless F is diagonal, 
which is rarely the case). 

^^ The term "window function" should not be confused with 
the term "EoR window". The former refers to the weights that 
specify the hnear combination of the true bandp owe rs that each 
bandpower estimate represents, as per Equation l |ll| . The latter 
refers to the region on the k^-k« plane that naturally has very low 
levels of foreground contamination, as illustrated in Figure 111 
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As a compromise option, we advise M ^ F^i/^ (again 
after a normalization of each row so that the window 
functions sum to unity). This choice gives window func- 
tions that are narrower than those for a diagonal M while 
maintaining reasonably small error bars. In addition, an 
inspection of Equation ^ reveals that this method gives 
a diagonal S, which means that errors between different 
bandpowers are uncorrelated. 

In Section [3^ we use MWA data to demonstrate the 
crucial role that the M ^ F~^/^ choice plays in preserv- 
ing the cleanliness of the EoR window. 

2.7. A real-world obstacle: ensuring that binning doesn't 
destroy error properties 

In previous sections, we have discussed how one can 
preserve all the desirable properties of the power spec- 
trum estimator of Section 12.11 in the face of all the real- 
world complications presented in Sections |2.2| through 



2.6 The result is an optimal yet practical estimator 
for the cylindrical power spectrum Pcyi{k±, /cy). We now 
turn to the problem of binning the cylindrical power spec- 
trum into the cosmologically relevant spherical power 
spectrum Psph(fc)j with a special emphasis on the preser- 
vation of the optimal properties of our estimator. 

Just as with the cylindrical power spectrum, we param- 
eterize the spherical power spectrum as piecewise con- 
stant, so that all the information is encoded in a vector 
of bandpowers p'^P^, so that: 



/cf' 



^sph (, fS 



(32) 



The spherical bandpowers are related to estimates of the 
cylindrical bandpowers p'^^' by the equation 



^cyl 



Ap' 



sph 



e, 



(33) 



where A is a matrix of size A'cyi x Nsph of Is and Os that 
relates fc_L-fc|| pairs to k bins, with A"cyi and A'sph equal to 
the number of cells in the fc_L-fc|| plane and the number of 
spherical k bins respectively. The vector £ is a random 
vector of errors on p'^''^ It has zero mean (assuming that 
one has taken the care to avoid additive bias in our esti- 
mator of the cylindrical bandpowers, as discussed above), 
but non-zero covariance equal to 'S'^^ = (ee*), where 
S'^-*' is given by either Equation (231 or (24), depending 



on whether the cylindrical bandpowers were computed 
using cross or auto-power spectra. (The methods pre- 
sented in this section are applicable either way). 

Our goal is to construct an optimal, unbiased estimator 
of p'^P^ from p'^y'. This is a solved problem (Tegmark 



1997a), and the best estimator p*'?'' is given by 



:>sph 



A*S-/iA]-iA*S-ir^' 



(34) 



with the final error covariance on the spherical bandpow- 
ers given by 



-isph 

-'aP 



(pL^'pJ') 



^^T){pt) 



A*E->]-i. (35) 



Since the A matrix has (by construction) a single 1 per 
row and zeros everywhere else, an inspection of Equa- 
tion ( 35 ) reveals that a diagonal S'^-^ implies a diagonal 
S^'P . In other words, the estimator given by Equation 
(34) preserves the decorrelated nature of the M ~ F^^/^ 



version of the cylin drical power spectrum estimator de- 
fined in Section lTBl This will not be the case for an arbi- 
trary estimator (such as one that is formed from taking 
uniformly weighted Fast Fourier Transforms, then squar- 
ing and binning) . We also emphasize that if one does not 
choose to u se d ecorre late d cylindrical bandpower vectors. 
Equations ( 34 ) and ( 35 1 require that one keep full track 



of the off-diagonal terms of S^ j. Without it, a consistent 
propagation of errors to the spherical power spectrum is 
not possible, and may even lead to a mistakenly claimed 
dete ction of the cosmological signal, as we discuss in Sec- 
tion |331 

Just as with the cylindrical power spectra, we would 
like to compute the window functions. The definition of 
the spherical window funct ions are exactly analogous to 
that provided in Equation (11 1 for the cylindrical power 
spectrum, so that 



^psph^ = W^P'^p'^P^. 



(36) 



Taking the expectation value of Equation ( 34 1 , we have 



(p^ph) = [A*S->]-iA*E-^(p=yi) 

= [A*S-ylA]-lA*E-lW=>'lAp^P^ (37) 

where we have used the definition of the cylindrical win- 
dow functions to say that (p'^^') = Wp^^^^ as well as the 
fact that p'^^' = Ap'^P'^ (with no error term because we 
are relating the true cylindrical bandpowers to the true 
spherical bandpowers). Inspecting this equation, we see 
that 

^sph ^ [A*-E^^\A]-^A*-E^^\W''y^A. (38) 

Therefore, by measuring the width of the spherical win- 
dow functions (rows of W'^p'^), one can place rigorous 
horizontal error bars on the final spherical power spec- 
trum estimate. 

2.8. Summary of the issues 

In the last few sections, we have provided techniques 
for dealing with a number of real-world obstacles. These 
include: 

1. Taking advantage of the flat-sky approximation 
and the rectilinearity of data cubes, as well as the 
conjugate gradient algorithm for matrix inversion 
to allow large data sets to be analyzed quickly. 

2. Using cross-power spectra rather than auto-power 
spectra in order to eliminate noise bias. 

3. Replacing inverses with pseudoinverses to deal with 
data that has missing spatial modes (due to incom- 
plete uv coverage) and missing frequency channels 
(due to RFI). 

4. Performing power spectrum decorrelation to avoid 
the leakage of foreground power into the EoR win- 
dow. 

5. Binning of cylindrical power spectra into spherical 
power spectra in a way that preserves desirable er- 
ror properties. 
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Crucial to this is the lact that these techniques aU oper- 
ate under a self-consistent framework. This allows faith- 
ful error propagation that accurately captures how var- 
ious real - world effects act t ogether. For example, it was 
shown in Dillon et al. ( 2013 1 that properly accounting for 
pixelization cttects in Equation ([5]) results in low Fisher 
information at high kn , providing a marker for parts of 
the fc_L-fc|| plane that cannot be well-constrained because 
of finite spectral resolution. The identification of such a 
region would be trivial if one had spectrally contiguous 
data, for then one would simply say that the largest mea- 
surable fell was roughly 1/ALu, where ALu is the width 
of a single frequency channel mapped into a cosmological 
line-of-sight distance. However, such a straightforward 
analysis no longer applies when there are RFl gaps in 
the data at arbitrary locations. In contrast, the unified 
framework presented in this paper allows all such com- 
plications to be folded in correctly. 

3. A WORKED EXAMPLE; EARLY MWA DATA 

Now that we have bridged the gap between theoreti- 
cal techniques for analyzing ideal data and the numer- 
ous challenges presented by real data, we are ready to 
bring together our methods, specify a covariance model, 
and estimate power spectra from MWA 32-tile prototype 
(MWA-32T) data. The data were taken between the 21st 
and 29th of March 2010, the first observing campaign 
during which data were taken that were scientifically use- 
ful. The observation s are described in more detail by 
Williams et al. (2012). Real data affords us two oppor- 
tunities. In this section, we look at the data to examine 
and quantify the differences between power spectrum es- 
timators and the pitfalls associated with choice of esti- 
mator. In Section HI we take advantage of everything we 
have developed to arrive at interesting new foreground 
results and a limit on the 21 cm brightness temperature 
power spectrum. 

3.1. Description of observations 

All of the data used for this paper were taken on 
the MWA-32T system. This system has since b een up 
graded to a 128-tile instrument (MWA-128T; [Tingay 
et al. ( 2013| )), but in this paper we focus exclusively on 
MWA-35T data, reserving the MWA-128T data for fu- 
ture work. 

The MWA-32T instrument consisted of 32 phased- 
array "antenna tiles" which served as the primary col- 
lecting elements. Each tile contained 16 dual linear- 
polarization wideband dipole antennas which were com- 
bined to form a steerable beam with a full width at half 
maximum (FWHM) size of ~ 25° at 150 MHz. The array 
had an approximately circular layout with a maximum 
baseline length of ~ 340 m, and a minimum baseline 
length of 6.6 m, although the shortest operating baseline 
during this observational campaign was 16 m. After dig- 
itization, filtering, and correlation, the final visibilities 
had a 1 second time resolution and 40 kHz spectral res- 
olution over a 30.72 MHz bandwidth. The instrumental 
capabilities are summarized in Table [l] 

For our worked example, we concentrate on March 
2010 observations of the MWA "EoR2" field. It 
is centered located at R.A.(J2000) = lO'' 20"" 0^ 
decl.(J2000) = -10° 0' 0", and is one of two fields at high 
Galactic latitude that have been identified by the MWA 



collaboration as candidates for deep integrations, ow 

ing to their low brightness temperat ure in low frequency 

measurements of Galactic em ission fHaslam et alj|1982 



de Oliveira-Costa et al.|[2008). For further details about 



the observational campaign or the EoR2 field, please see 
Williams et al.^ (2012]) , which was based on the same set 
of observations as the ones used in this paper. 

Observations covered three 30.72 MHz wide bands, 
centered at 123.52 MHz, 154.24 MHz and 184.96 MHz, 
corresponding to a redshift range of 6.1 < z < 12.1 
(the redshift range of the results presented in this work 
is slightly smaller because of data flagging) for the 
21 cm signal. The 123.52 MHz and 154.24 MHz bands 
were observed for approximately 5 hours each, and the 
184.96 MHz band was observed for approximately 12 
hours. 

These early data from the prototype have provided us 
with a set of test data that enabled development of ex- 
tensive analysis methods and software on which the re- 
sults of this paper are based. The early prototype had 
shortcomings (e.g., mismatched cables, receiver firmware 
errors, correlator timing errors) that compromised the 
calibration to some extent, raising the apparent noise 
level. Additionally, the instrument was only operating 
with < 29 tiles, and with a 50% duty cycle throughout 
the course of these observations. We account for this 
in Section [3. 3| by determining the magnitude of the noise 
empirically, in order to be able to place rigorously conser- 
vative upper limits on the cosmological power spectrum. 
We expect that data from later prototype campaigns and 
from the full array will produce result closer to theoreti- 
cal expectations. 

3.2. Mapmaking 

Before the data can be used as a worked example for 
our power spectrum estimator, however, we must convert 
the measured visibilities into a data cube of sky images 
at every frequency in our band. In other words, we must 
form the data vector x, defined by Equation (pi), which 
serves as the input for our power spectrum pipeline. 

To form the data vector, we performed the following 
steps. First, we perfo rmed a reductio n proce dure similar 



to that described in 
tial flagging and call 



(2012 



I 



for the ini- 
ydra A was 



Williams et al. 

Dration of the data 
identified as the dominant bright source in the field, and 
used for calibration assuming a point source model. The 
Hydra A source model was then subtracted from the uv 
data. As this same source model was also used for gain 
and phase calibration, this can be thought of as a "peel- 
ing" source removal procedure (iN oordam 2004; v an der 
Tol et aL|2007{|MitcheU et al.|20(te||Intema et al 1200 9^) on 



a single source. Alternatively, in the absence of gridding 
artifacts, this is equivalent to imaging the point-source 
model and subtracting it from the data as part of the 
direct foregrou nd subtraction step d iscussed in the first 
step of Section [23| ( |Tegmarkp997a[ ). 

The subtracted data were imaged using the CAS A 
task clean without deconvolution to produce "dirty" 
images. No multi-frequency synthesis was performed, 
so that the full 40 kHz spectral resolution of the data 
would be available. The visibilities were gridded using 
w-projection kernels ( Cornwell et al.|[2008 ) with natural 
(inverse- variance) weighting to produce maps at each fre- 
quency with a cell size of 3' over a 25.6° field of view. The 
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Table 1 

MWA-32 Instrument Parameters 

Field of View (Primary Beam Width) ~ 25° at 150 MHz 

Angular Resolution 

Collecting Area 

Polarization 

Frequency Range 

Instantaneous Bandwidth 

Spectral Resolution 



' 20' at 150 MHz 
~ 690 m^ towards zenith at 150 MHz 
Linear X-Y 
80 MHz to 300 MHz 
30.72 MHz 
40 kHz 



resulting cubes contained ~ 200 million voxels, with 512 
elements along each spatial dimension and 768 elements 
in the frequency domain. It is important to note that the 
pre-flagging performed on the data resulted in the flag- 
ging of entire frequency bands (which means that there 
are gaps in the final data cube). Cubes were generated 
for each 5 minute snapshot image. 

The individual snapshot data cubes were combined us- 
ing the prim ary beam inverse- varian ce weighting method 



described in Williams et al. (^2012). The weighting and 



primary beams were simulated separately for each 40 kHz 
frequency channel in each 5 minute snapshot. The com- 
bined maps and weights were saved, along with the effec- 
tive point spread function at the center of the field. Two 
additional data cubes were created by averaging alternat- 
ing 5 minute snapshots (i.e. even numbered snapshots 
were averaged into one cube, and odd numbered snap- 
shots were averaged into the other) so that they were 
generated from independent data, but with essentially 
the same sky and uv coverage properties. 

A further flux scale calibration of the integrated 
cubes was performed using three bright point sources: 
MRC 1002-215, PG 1048-090, and PKS 1028-09 to set 
the flux scale on a channel-by-channel basis. A two di- 
mensional Gaussian fitting procedure was used to fit the 
peak flux of each of these sources in each 40 kHz chan- 
nel of the data cube. Predictions for each source were 
derived by fltting a power law to source measur ements 
from the 4 85 GHz Parkes-MIT-NRAO survey (Grifiith 
et al.||1995|) , the 408 MHz Molongio Reference C atalog 
IfCargeeTa l. 1981), the 365 MHz Texas Survey (|Dou 



glas et aL||1996|), the 160 MHz and 80 MHz Culgoora 



Source List (|Slee||1995| and the 74 MHz VLA Low 



frequency Sky Survey (Cohen et al. 20071. A weighted 



least-squares fit was then performed to calculate and ap- 
ply a frequency-dependent flux scaling for the cube to 
minimize the square deviations of the source measure- 
ments from the power law models. 

An additional flagging of spectral channels was per- 
formed based on the root-mean-square (RMS) noise in 
each spectral channel of the cube. A smooth noise model 
was determined by median flltering the RMS channel 
noise as a function of frequency (bins of 16 channels were 
used in the filtering). Any channel with 5cr or larger devi- 
ations from the smoothed noise model was flagged. Upon 
inspection, these additional flagged channels were ob- 
served to be primarily located at the edges of the coarse 
digital filterbank channels, which were corrupted due to 
an error in the receiver firmware. After this procedure, 
approximately one third of the spectral channels were 
found to have been flagged. 

Each individual map covered 25.6° x 25.6° at a res- 
olution of 3' with 768 frequency channels (40 kHz fre- 



quency resolution). To decrease the computational bur- 
den of the covariance estimation, each map was subdi- 
vided into 9 subflelds, and the pixels were averaged to 
a size of 15'. The data cubes were mapped to comov- 
ing cosmological coordinates using WMAP-7 derived cos- 
mological parameters, with J^m — 0.26 6, SIa = 0.734, 



Hn = 71 km s Mpc , and fl^ = (Komatsu et al 



2011| ). 

At this point, the data cubes were ready to be used as 
input data to our power spectrum estimator, i.e., we had 
arrived at the final form of the data vector x. However, 
estimating power spectra and error statistics using the 
formalism of Section |2] also requires a covariance model, 
which we construct in the next section. 



3.3. Covariance model 



We follow Liu & Tegmark (2011) and Dillon et al. 



( 2013 1 in modeling the covariance matrix C as the sum of 
independent parts attributable to noise and foregrounds. 
We leave off the signal covariance because it only con- 
tributes to the final error bars by accounting for cosmic 
variance — a completely negligible effect in comparison to 
foreground and noise-induced errors. We adopt a conser- 
vative model of the extragalactic foregrounds by treating 
them as a Poisson random field of sources with fluxes less 
than 100 Jy, after the manual removal of Hydra A. By 
treating all extragalactic foregrounds as "unresolved," we 
effectively throw out information about which lines of 
sight are most con taminated by bright foregrounds. As 
Dillon et al. (20131 showed, future analyses can improve 
on our limits by including more information about the 
foregrounds. We begin with the pa rameterized covari- 
ance model of Liu & Tegmark (2011 ), 



C" 



olvcd 



== 1.4 X 10 
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exp 
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-1 




1 ^ ^pix 

J V / 















(39) 



where Vi, = 150 MHz is a reference frequency, ui is the 
frequency of the ith voxel, which has an angular distance 
of r±i from the field center. The spectral index is R = 
0.5, the uncertainty in the spectral index is CTk = 0.5, 
the clustering correlation length is <t± = 7', flpi^ is the 
angular size of each pixel, the fiux cut S'cut = 100 Jy, and 
dn/dS is the differential source count from Di Matteo 
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eral] ( |2002| , 
dn 



dS 



(4000 Jy^^sr-i) 




-2.51 



-1.75 



for S > 0. 



for S < 0.. 



Jy 
Jy. 



(40) 



We adapt this model for the fast pow er spectrum esti- 
mation method outhned in Section [2?2] by calculating the 
translationally invariant approximation to t his model in 



the manner described in Dillon et al 
For the Galac tic synchr otron, 



& 



( |2013]) . 

also follow [Liu 
Tegmark] (I2011J) and Dillon et al.| ( |2013[ ) for the pa- 
rameterization of the synchrotron emission covariance. 
Namely, we adopt k = 0.8, a^, = 0.1, a± = 30°, and re- 
place the first three terms of the covariance in Equation 



39 



with T2y„,h 



(335.4 K)2. 

Our model for th e instrumental noise is adopted from 
Dillon et al. ( 2013[ ), with one key difference: the overall 
normalization. I'or each subband, we let the noise covari- 
ance matrix scale by a free multiplicative constant. This 
is equivalent to treating the combination Tgyg/(A^„^iobs) 
as a free parameter. We then fit for that parameter by re- 
quiring the RMS difference between the two time slices — 
which should be free of sky signal — for the densely sam- 
pled inner region of uv space and rescaling our noise co- 
variance matrix to match. The spatial structure of the 
covariance was left unchanged. Even though the data is 
somewhat nosier than suggested by a first principles cal- 
culation assuming fiducial values for system temperature 
and antenna effective area, this empirical renormaliza- 
tion allows for an honest account of the errors introduced 
by instrumental effects. 

To verify that our parameterization of the foregrounds 
was reasonable, we varied these parameters over an or- 
der of magnitude and found that they had little effect on 
our final power spectrum estimates, except at the lowest 
values of k. There are two reasons for this: first, since 
we are only measuring the power spectrum of the sky, we 
need not worry about precisely subtracting foregrounds. 
Second, because the noise in our instrument is still more 
than two orders of magnitude from the cosmological sig- 
nal, in the EoR window our band power measurements 
will be noise dominated and agnostic to our foreground 
model. Future analyses might include a more thorough 
treatment of the foregr ounds, especially by utilizing the 



full power of the Dillon et al. (2013) method to include 



information about the positions, ffuxes, and spectral in- 
dices of individual point sources. 

3.4. Evaluating power spectrum estimator choices 

With both a data vector x and a covariance matrix 
C in hand, we can now apply the methods of Section 
[2] to estimate power spectra. In doing so, we deal with 
real-world obstacles using all of the techniques that we 
have developed. In this section, we show why all this is 
necessary. 

In Section [2. 6| we touted the choice of power spectrum 
estimator p = Mq with M ~ F~^/^ as a compromise 
solution between the choice with the smallest error bars, 
M ~ I, and the choice with the narrowest window func- 
tions, M ~ F^^. In the race to detect the power spec- 



trum from the EoR, one might be tempted to aggressively 
seek out the smallest possible errors. This could prove a 
deleterious choice, as we will now show using MWA-32T 
data. 

First, in Figure [2] we compare cylindrical power spec- 
tra, p, generated using two different estimat ors of the 

■ ^ On 
with 



in^ 



power spectrum that we presented in Section 
the left, we have used M ~ I, the estimator 
the smallest error bars, and on the right we have used 
M ~ F"-'^/^, the estimator with decorrelated errors. In 
both cases, we have plotted the absolute value of the 
power spectrum estimates (which can be negative be- 
cause they are cross-power spectra). Because the two 
estimates are related to one another by an invertible ma- 
trix, they contain the same cosmological information. In 
a sense, the M ~ F~^/^ method is the most honest es- 
timator of the power spectrum because the band powers 
form a mutually exclusive and collectively exhaustive set 
of measurements. In other words, they represent all the 
all the power spectrum information from the data, di- 
vided into independent pieces. 

Moreover, just because two sets of estimators have the 
same information content does not mean that they are 
equally useful for distinguishing the cosmological power 
spectrum from foreground contamination. In Figure [2J 
the minimum variance estimator for the power spectrum 
introduces considerable foreground contamination into 
the EoR window, demarcated by the expected angular 
extent of the wedge feature (which we introduced in Sec- 
tion IT] and will discuss in greater detail in Section 4.1 1. 
Eveiinighly suspect features at high k± where uv cover- 
age is spottiest seem to get smeared across k± and into 
the EoR window. We cannot simply cut out the wedge 
from our cylindrical-to-spherical binning and expect a 
clean measurement of the power spectrum in the EoR 
window. 

Looking closely at Figure[2J one might notice that some 
regions of the EoR window on the lefthand panel still 
seem very clean — cleaner perhaps that the same regions 
in the righthand panel. To examine that apparent fact, 
we plot Pa instead of \pa \ in Figure pj To make the 
figure more intelligible, we have plottea colors based on 
an sinh^^ color scale with a sharp color division at 0. 
The sinh"^ has the advantage of behaving linearly at 
small values of p^ and logarithmically at large positive 
or negative pa . 

What emerges is a striking difference between the two 
For the reasons discussed in Section 12.31 we 



estimators. 

have chosen to estimate the cross power spectrum be- 
tween two time-interleaved sets of observations. As a re- 
sult, we expect that instrumental noise should be equally 
likely to contribute positive power as it is to contribute 
negative power. In noise dominated regions of the k±-k^\ 
plane, we expect about half of our measurements to be 
positive and about half to be negative. That is exactly 
what we see in the EoR window of the M ^ F~^/^ es- 
timator. However, the M ~ I estimator in the lefthand 

^^ In our comparison of choices for M, we drop the M ~ F~^, 
5-function windows choice. In addition to proving the noisiest es- 
timator, it suffers from strong anti-correlated errors. We adopt 
the perspective that the important comparison is between the "ob- 
vious" choice, the minimum variance M ~ I, and our preferred 
choice with decorrelated errors, M ^ F^^' ^. 
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Figure 2. Unless one chooses a power spectrum estimator with decorrelated errors, foregrounds and other instrumental effects can leak 
significantly into the EoR window. Here we show the absolute value of the cylindrical power spectrum estimate from the subband centered 
on 158 MHz {z = 8.0) and averaged over all 9 fields. On the left, we have set M ^ I. On the right, M ~ F~^' ^. We expect contamination 
from smooth spectrum foregrounds interacting with the chromatic synthesized beam to occupy the "wedge" portion of Fourier space, defined 
in Equation ([l]l. Optimistically, the wedge is delimited by the extent of the main lobe of the primary beam; conservatively, we should not 
see bright foreground contamination beyond the horizon. In the regions where the power spectrum is noise dominated, we expect little 
structure in the fcii direction in the EoR window above some moderate value of fcii. In the left panel, we see considerable ka structure in 
the form of horizontal bands, attributable to foreground contamination and instrumental effects, that has leaked into the putative EoR 
window. 



panel clearly shows positive power throughout the entire 
supposed EoR window. Though the magnitude of that 
power is not enormous — often it is well within the ver- 
tical error bars — the overall bias towards positive cross 
power means that sky signal is contaminating the EoR 
window. This is pre cisely the problem we were worried 
about in Section [2.6| and the data have clearly manifested 
it. 

This also explains why there appeared to be less power 
in the EoR window of the lefthand panel of Figure 2j by 
taking the absolute value of the weighted average ofpos- 
itive and negative quantities, we expect to measure a 
smaller absolute value of the power. However, as this 
figure clearly shows, that weighted average is biased by 
foreground leakage. And, even though there still appears 
to be a region just inside the EoR window that retains 
positive band powers consistent with foregrounds, that 



small amount of leakage can be attributed to finite sized 
windows functions and to calibration uncertainties. Re- 
gardless, it does not appear to be an insurmountable lim- 
itation to the cleanliness of the EoR window; rather, it 
suggests that we should be careful in how we demarcate 
the EoR window when calculating spherically-averaged 
power spectra. 

In addition to producing a cleaner EoR window, the 
decorrelated estimator of the power spectrum yields an- 
other advantage: narrower window functions. Both the 
estimator with the minimum variance and estimator with 
decorrelated errors represent, in aggregate, the weighted 
average of the true, underl ying band power spectrum. 

In Figure |4] we show 



2.1 



as we discussed in Section 

the improvement that the decorrelated estimator offers 

over the minimum variance estimator by narrowing the 
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Figure 3. One advantage of caleulating the cross power spectrum of interleaved time-slices of data is that we can easily tell which regions 
of Fourier space are noise dominated. Here we reproduce the power spectra from Figure [2] without taking the absolute value of P{k). By 
plotting with a discontinuous, sinh"'^ color scale, it is easy to see that the EoR window for our decorrelated power spectrum estimate (right 
panel) has roughly an equal number of positive and negative band power estimates — exactly what we would expect from a noise dominated 
region. By contrast, our power spectrum estimate with correlated errors (left panel) shows positive power over almost all of Fourier space, 
indicating ubiquitous leakage of contaminants into the EoR window. 



window functions considerablyPj We show five example 
window functions from the same subband that we plot in 
Figure [2J cropped to fit together on one set of axes, each 
centered at their respective peaks. Because the window 
functions are normalized to sum to 1, the breadth of each 
window function is reflected by the value of the central 
peak. As we expected, the window functions are con- 
siderably narrower for our decorrelated power spectrum 
estimator. 

Even after binning from cylindrical power spectra to 
spherical power spectra, the difference remains quite 
stark. In Figure [5 we see clearly that choosing a power 
spectrum estimator with decorrelated errors also consid- 
erably improves the window functions in one dimension 

^^ While the choice of M ^ F~^ " ensures that the power spec- 
trum estimator covariance is diagonal (recall, S = MFM* while 
W = MF), it does not mean that the window functions are delta 
functions. The off-diagonal terms of S might be zero even if the 
off-diagonal terms of W are not. 



as well as two. 



In Section 2.7 we argued that an inverse covariance 
weighted binning scheme for estimating spherical band 
powers produced optimal spherical power spectrum esti- 
mate. In the case where M is chosen either for the small- 
est possible error bars or the narrowest possible window 
functions, the estimator covariance S'^^ is non-diagonal. 
Assuming that the matrix is actually diagonal, at one or 
more steps in the binning and error propagation, can lead 
error bars that are overly conservative or — worse yet — 
error bars that are insufficiently conservative and might 
falsely lead to a claimed detection. In Figure |6j we show 
the effects of making a suboptimal choice for binning. 

If one fully models the covariance matrix Yf^ , includ- 
ing off-diagonal elements, but chooses to generate ^^^^ as 
an inverse variance (and not inverse covariance) weighted 
average of cylindrical band powers, neglecting off diago- 
nal terms in the weighting, one's estimators will be nois- 
ier as a result (see the solid lines in Figure [6]) . These are 
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Figure 4. By using an estimator of the 21 cm power spectrum with uncorrelated errors, we significantly narrow the window functions 
that relate the ensemble average of our estimator to the true, underlying power spectrum. Here we show a sample of five cropped window 
functions for the power spectrum estimate in Figure [2] each centered at their maxima, for both an estimator with correlated errors (left 
panel) and an estimator with uncorrelated errors (right panel). Though the estimator with correlated errors produces smaller vertical 
error bars, it acheives this by "over-smoothing" many band powers together. Narrow window functions let us independently measure many 
modes of the power spectrum. The band power measured with M ~ F~^''^ is one of a set of mutually exclusive and collectively exhaustive 
pieces of information. 



the correct errors for the suboptimal choice of estimators. 

Even worse, if one assumes that 11'^^ is diagonal when 
it is not, one is led either to overestimate the error bars, 
in the case of M '-^ F~^, or underestimate them, as would 
be the case when M ~ I. This is because the former 
case general exhibits anti-correlated errors while the lat- 
ter suffers from correlated errors. The last scenario is the 
most troubling: by aggressively choosing the estimator 
with the smallest vertical error bars (M ^ I) and then 
neglecting the correlations between errors, one will un- 
derestimate the error bars and might be lead to falsely 
claiming a detection. In this case, the estimator is sub- 
optimal and the errors are incorrect. Additionally, as we 
show in Figure [Tj if one were to calculate the the window 
functions under the assumption that S*^^ is diagonal, 
one would find window functions several times boarder 
than they would otherwise be. 

Thankfully, choosing the cylindrical power spectrum 



estimator with decorrelated errors avoids the subtle dif- 
ference between inverse variance and inverse covariance 
weighted binning. The M ~ F~^/^ decorrelated estima- 
tor preserves the EoR window and allows for easy, op- 
timal binning of uncontaminated regions into spherical 
band power spectrum estimates. 

4. EARLY RESULTS 

Having developed and demonstrated a technique that 
robustly preserves the EoR window while thoroughly and 
honestly keeping track of the errors on and correlations 
between our band power estimates, we can now confi- 
dently generate some interesting preliminary science re- 
sults. Because these data span the widest redshift range 
to date, we are able to investigate the behavior of the 
wedge feature over many frequencies. Understanding the 
behavior of the EoR window over a large redshift range 
is important, since there is still considerable uncertainty 
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Figure 5. Even after optimally binning the cylindrical power 
spectra from Figure p] to spherical power spectra, the choice of a 
power spectrum estimator with decorrelated errors produces much 
narrower window functions than the minimum variance technque. 
In addition to maintaining a clean EoR window, the choice of 
M ~ F~^ " provides the additional benefit of allowing power spec- 
trum modes to be measured more independently. 



about the timing and duration of the EoR. Moreover, 
it is often argued that a tentative first detection of the 
cosmological signal wih only be convincingly distinguish- 
able from residual foregrounds if one can show that the 
21cm brightness temperature fluctuations peak at some 
redshift, since theory predicts that the midpo int of reion- 
ization should be marked b y such a peak (Lidz et al. 
Bittner & Loeb 20111). It is therefore essential to 



foregrounds) over a broad frequency range. We also ap- 
ply our methods from Section [2] to calculate spherically 
averaged power spectra over our entire redshift range, 
including error bars and window functions, thus setting 
a limit on the 21 cm brightness temperature power spec- 
trum during the EoR. 

4.1. The wedge 

In Figure [H we show all the cylindrical power spectra 
over the redshift range probed by our current observa- 
tions. The spectra are sorted into three rows, each of 
which contain data coming from a single 30.72 MHz wide 
frequency band. All of the spectra were generated using 
the same techniques that were used to generate the ex- 
ample cylindrical power spectra in Section |3.4| and thus 
contain all the desirable statistical properties discussed 
in Section |3.4| One sees that in every case the fore- 
grounds are mostly confined to the wedge region in the 
bottom right corner of the fc_L-fc|| plane. This builds upon 



2008 



characterize the EoR window (and by extension, residual 



the single frequency observations of Pober et al. (2013), 
demonstrating the existence of the EoR window across a 
wide range of frequencies relevant to EoR observations. 

Having these measurements also allows us to exam- 
ine the behavior of the EoR window as a function of 
frequency. Consider first the high k± regions of the k±- 
fcjj plane. The most striking feature here is the wedge. 
Consistent with being dominated by foreground power, 
the wedge generally gets brighter with decreasing fre- 
quency within each wide frequency band, just as fore- 
ground emission is known to behave. The extent of the 
wedge is also in line with theoretical expectations. Re- 
call from Equation (llj that the wider the field-of-view, 
the farther up in fcy the wedge goes. Since the field-of- 
view is defined by the primary beam, whose extent de- 
creases with increasing frequency, one expects the wedge 
to have the largest area at the lowest frequencies. This 
trend is clearly visible in the cylindrical power spectra 
of Figure Isj where the wedge extends to the highest fcy 
at the hignest redshifts. Importantly, the wedge is con- 
fined to its expected location across the entire range of 
the observations. To see this, note that we have overlaid 
Equation (fTl) on the plots, with the dashed line corre- 
sponding to ^max equal to that of the first null of the 
primary beam, and the solid line with ^max = '"'/Z (the 
horizon). At all frequencies, the most serious contami- 
nations lie within the first null, ensuring that the EoR 
window is foreground-free. 

Foregrounds also enter indirectly into the instrumental 
noise-dominated regions because the MWA is sky-noise 
dominated. Thus, as the brightest sources of emission in 
our observations, the foregrounds set the system temper- 
ature, and result in a higher instrumental noise at higher 
redshifts. This trend can be seen within each wide fre- 
quency band (each row of Figure Isl), although the shght 
interruption of this trend between oands suggests an ad- 
ditional source of noise. 

At low k±, theory suggests that foregrounds will con- 
taminate a horizontally-oriented region at the bottom of 
the plot. This is clearly seen in the highest frequency 
plots. Interestingly, at lower frequencies the increasing 
instrumental noise plays more of a role, and the fore- 
ground contribution is less obvious in comparison (al- 
though it is still there). While a naive reading of some 
of these low frequency plots (such as the one for z = 9.1) 
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Figure 6. Neglecting the fact ttiat the covariance of the power spectrum estimator is, in general, non-diagonal, can lead to two mistakes 
that can either unnecessarily enlarge our error bars or, even worse, unjustifiably shrink them. In this figure, we first show an approximately 
10% increase in the vertical error bars on the power spectrum ( soli d lines) from a suboptimal inverse variance weighted binning scheme, 
rather than the inverse covariance weighted binning of Equation | |34[ l . This problem is obviated by choosing an estimator with decorrelated 
errors and thus a diagonal covariance matrix. If one simply assumes that the estimator covariance in Equation ( |35| is diagonal when it is 
not (dotted lines), one is led, depending on the choice of estimator, either to roughly 50% larger error bars than necessary or, worse yet, 
artificially small error bars. The last mistake, choosing an estimator with small error bars — despite its wide window functions — and then 
neglecting the off-diagonal terms in the estimator covariance, is potentially the most pernicious since it could lead to a claimed detection 
in the absence of signal. 



might suggest that the EoR window extends to the low- 
est fc||, such a concl usion would be misguided. As we 
shall see in Section 4.2[ these modes are likely dom- 
inated by foregrounds (and therefore do not integrate 
down with further integration unlike instrumental noise 
dominated modes). Moreover, the error statistics (which 
self-consistently include foreground errors in our formal- 
ism) suggest that low fcy modes are less useful for con- 
straining theoretical models, and that the true EoR win- 
dow does in fact lie at higher fcy, as suggested by the- 
ory. Again, this highlights the importance of estimating 
power spectra in a framework that naturally contains a 
rigorous calculation of the errors involved. 



4.2. Spherical Power Spectrum Limits 

Having confirmed that the EoR window behaves as ex- 
pected, we will now proceed to place constraints on the 
spherical power spectrum. In top panel of Figure [9] we 
show the result of binning the z = 10.3 cylindrical power 
spectrum of Figure [Sl usin g the optimal binning formulae 
presented in Section |2.7| In addition, for ease of inter- 



pretation, we elect to plot 



A(fc) 



2n 



,P{k) 



(41) 



(which simply has units of temperature) rather than 
P{k) itself. 

To quantify the errors in our spherical power spec- 
trum estimate, we also bin the cylindrical power spec- 
trum measurement covarianc es a nd window functions us- 
ing the formulae of Section 2.7 The resulting window 
functions are shown in the bottom panel, and give an 
estimate of the horizontal error bars. Thinking of these 
window functions (which, recall, are normalized to inte- 
grate to unity) as probability distributions, the horizon- 
tal error bars shown in the top panel are demarcated by 
the 20th and 8Gth percentiles of the distribution. (This 
corresponds to the full-width-half-maximum in the event 
that the window functions are Gaussians). The vertical 
error bars were obtained by taking the square root of 
each diagonal element of the covariance matrix. Since 
the methods of Section [277| carefully preserved the diago- 
nal nature of the bandpower covariance, each data point 
in Figure |9] represents a statistically independent mea- 
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Figure 7. Just as with the error bars in Figure [6] generating 
suboptimally binned spherical power spectrum estimates by ne- 
glecting off-diagonal terms in the estimator covariance can lead 
to wider window functions than necessary. We illustrate the ef- 
fect by comparing the width of the window functions between the 
20th and 80th percentiles between the two binning schemes. This 
is important for the choice of power spectrum estimator with the 
smallest error bars and widest window functions (M ~ I). In the 
case where our power spectrum estimator has uncorrelated errors, 
there are no off-diagonal terms in the estimator covariance and 
both binning schemes are identical. In the case of the estimator 
with <5-function window functions, suboptimal binning does not af- 
fect the window functions — though it still affects the vertical errors 
(see Figure mjl. 

surement. This would not have been the case had we 



not employed the decorrelation technique of Section 2.6 



Immediately obvious from the plot is that there is 
qualitative difference between the data points at low k 
and those at high k. In particular, the points at low k 
are detections of the sky power spectrum, whereas the 
points at high k are formal upper limits. This is not to 
say, of course, that the cosmological EoR signal has been 
detected at low k. Rather, recall from Section |2.3| that 



in an attempt to avoid having to make large bias sub- 
tractions, we elected to compute cross-power spectra of 
total sky emission rather than of the cosmological si gnal , 



with the expectation (largely confirmed in Section |4.1 1 
that the intrinsic cleanliness of the EoR window would be 
sufficient to ensure a relatively foreground-free measure- 
ment at high fc|| . Now, our survey volume is such that 
we are sensitive almost exclusively to regions in Fourier 
space where fcy ^ k±. When binning along contours of 
constant k in the cylindrical Fourier space, we have that 

k 



k\ 



kl 



II , and therefore the low k points of 

Figure[9|map to low fc|| . The detections seen at low k thus 
reside outside the EoR window and are almost certainly 
detections of the foreground power spectrum. 

Despite the fact that the low k modes are foreground 
dominated, they still constitute a formal upper limit on 
the cosmological power spectrum, since the foreground 
power spectrum is necessarily positive. In fact, our cur- 



rent, most competitive upper limit resides at the lowest 
k values. However, this is unlikely to continue to be 
the case as more data is taken with the MWA, for two 
reasons. First, as foreground-limited measurements, the 
data points at low k will not average down with fur- 
ther integration time. In addition, the error statistics in 
the region are not particularly encouraging. The window 
functions (and therefore the horizontal error bars) are 
seen to broaden towards lower fc, reducing the ability of 
constraints at those k to place limits on theoretical mod- 
els. (This is most easily seen by recalling that the window 
functions integrate to unity by construction, and thus the 
increase in their peak values towards higher k implies a 
broadening of the window functions) . The broadening of 
the window function s is an ex pected consequence of fore- 
ground subtraction (Liu & Tegmark 2011) and thus will 
likely continue to liniit the usefulness of the low k regime 
unless future measurements can characterize foreground 
properties with exquisite precision. 

In contrast, the points at high k do reside in the EoR 
window. The constraints here are limited by thermal 
noise, as we saw in Section |3.4[ Bolstering this view is 
the fact that the data here are consistent with zero, as 
one expects for a noise-dominated cross-power spectrum. 
The limits her e ar e given by the 2cr errors pr edic ted by 
the Equation ( 35 1 . 

larger 



As mentioned in Section 3.3 



these 
e pre- 



errors are somewhat larger than what might 
dieted by a theoretical sensitivity calculation. However, 
they are consistent with rough estimates of the errors 
obtained from a calculation of root-mean-square values 
from the images produced in Section |3.2[ This suggests 
that the larger-than-expected errors are due to noisier- 
than-expected input maps, and not to any approxima- 
tions made in the power spectrum estimation techniques 
presented in this paper. The data on which these re- 
sults are based are from the very first operation of the 
prototype array, and we expect better performance in 
later data. Encouragingly, we note also that as noise- 
dominated constraints, the measurements at high k will 
continue to imp rove with integration time. 

In Figure |10[ we show power spectrum limits across 
the entire frequency range of the MWA, along with some 
theoretica.1 pre dictions generated using the models in 



Barkana (20091. At the lowest redshift, no theory curve 
is plotted because the model predicts that reionization 
is complete by then. This yet again underscores the 
importance of making measurements over a broad fre- 
quency range — with access to a wide range of redshifts, 
future detections of the cosmological signal can be dis- 
tinguished from residual foregrounds by measuring null 
signals at redshifts where rei oniz ation is complete. 

Each redshift bin of Figure 10 exhibits trends that are 
qualitatively similar to those discussed above for the z = 
10.3 case. Over all bands, our best limit is A(/c) < 0.3 K, 
occurring at z — 9.5 and k — 0.046 cMpc^^. However, as 
remarked in Section [3. 3[ the lowest k bins can be rather 
sensitive to the covariance model, and if one excludes 
those bins, our best limit is A(fc) < 2K, at z = 9.5 and 
k = 0.134cMpc~"'^. While our limits may no t be quite as 
low as other existing limit s in the literature ( Paciga et al. 
2013 Parsons et alj2013), they are the only limits on the 



EoR power spectrum that span a broad redshift range 
from z = 6.2 to z = 11.7. Moreover, these statistically 
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Figure 8. Examining cylindrically binned power spectra for each subband (each averaged over all nine subfields), reveals several important 
trends with frequency of the EoR window and the foregrounds. Each row is a single simultaneously observed frequency band. Since different 
bands were observed for different amounts of time, direct comparisons between rows is challenging. However, several clear trends emerge. 
For each band, moving to higher redshift (increasing wavelength) shows stronger foregrounds, a larger wedge (in part due to a wider 
primary beam), and a noisier EoR window (due to a higher system temperature). In general the brightest foreground contamination is well 
demarcated by the wedge line in Equation ^\ for the primary beam (dotted line) and especially by the wedge line for the horizon (solid 
line). In short, the wedge displays the theoretically expected frequency dependence. 
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Figure 9. Our method allows for the estimation of the spherically binned power spectrum in temperature units, A(fc), while keeping full 
acount of both vertical error bars and window functions (horizontal error bars) and making an optimal choice in the tradeoff between the 
two. In the top panel, we have plotted our spherical power spectrum estimates of the subband centered on 158 MHz (z = 8.0), including 
1(T errors on detections (which are often only barely visible), 2cr upper limits on non-detections, and horizontal error bars that span the 
middle three q uintiles of the window functions (bottom panel). At low fc, the wide error bars are the expected consequence of foreground 
contamination ( |Liu fc Tegmark|2011| l. Downward arrows represent measurements consistent with noise at the 2a level. Even though the 
area under the primary beam wedge tias been excised from the 2D-to -lD binning, the d etection of foregrounds at low k, is expected due to 
the contribution of unresolved foregrounds over a wide range of k± ( [Dillon et al.|2013^ . Our fiducial theoretical power spectrum is taken 
from [Barkana, (2009) . 
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rigorous limits will likely improve with newer and more 
sensitive data from the MWA. 

5. CONCLUSIONS 

In this paper, we have accomplished three goals. First, 
we adap ted 21 cm pow er spectru m estimation techn iques 
from Liu fc Tegmark| (J2011J) and [Dillon et al!] ( |2013[ ) with 
real- world obstacles in mmd, so that they could be ap- 
plied to real data. With early MWA data, our gener- 
alized formalism was then used to demonstrate the im- 
portance of employing a statistically rigorous framework 
for power spectrum estimation, lest one corrupt the nat- 
urally foreground-free region of Fourier space known as 
the EoR window. Finally, we used the MWA data to set 
limits on the EoR power spectrum. 

In confronting real- world obstacles, our desire is to pre- 
serve the statistical rigor in previous matrix-based power 
spectrum estimation frameworks. To avoid having to 
perform direct subtractions of instrumental noise biases, 
we advocate computing cross-power spectra between sta- 
tistically identical subsets of the data (in the case of the 
MWA worked example of this paper, these subsets were 
formed from odd and even time samples of the data). 
This has the effect of eliminating noise bias in the power 
spectrum, although instrumental noise continues to con- 
tribute to the error bars. To avoid direct subtractions 
of foreground biases, we simply look preferentially in the 
EoR window, where foregrounds are expected to be low. 
Missing data, whether from incomplete uv coverage or 
RFI flagging, can be dealt with using the pseudoinverse 
formalism. Doing this allows the effects of missing data 
to be self-consistently propagated into error statistics 
such as the power spectrum covariance and the window 
functions. In an effort to preserve the cleanliness of the 
EoR window, one should form decorrelated bandpower 
estimates, which have uncorrelated errors and reason- 
ably narrow window functions. Care must then be taken 
to preserve these nice properties via an optimal binning 
of cylindrical bandpowers into spherical bandpowers. 

Using early MWA data to demonstrate these tech- 
niques, we have confirmed theoretical predictions for the 
existence of the EoR window and have extended previous 
observations done by other groups to a much wider fre- 
quency range. This allowed us to check predicted trends 
of the EoR window as a function of frequency, all of 
which are consistent with theory. Crucially, we found 
that wit hout using the decorrelation technology of Sec- 
tion 2.6 measurements in the EoR window are not in fact 



instrumental noise dominated, and contain a systematic 
bias that is indicative of foreground leakage from outside 
the EoR window. 

The early MWA data has also allowed us to place lim- 
its on the cosmological EoR power spectrum. Our best 
limit is A(fc) < 0.3 K, at ^ = 9.5 and k = 0.046 cMpc 
(or A{k) < 2K at z = 9.5 and fc = 0.134 cMpc"^ if one 
discards the lowest k bin to immunize oneself against 
foreground modeling uncertainties). This may not be 
competitive with other published observations, but gen- 
eralizes them in an important way; instead of focusing 
on one particular frequency, our limits span a wide range 
of redshifts relevant to the EoR, going from z = 6.2 to 
z — 11.7. In addition, these limits will almost certainly 
improve in the near future, using already-collected (but 
yet to be analyzed) data from the MWA-32T system, as 



well as soon-to-be-collected data from the MWA-128T 
system. The rigorous statistical tools developed in this 
paper should be equally applicable to these newer data 
sets, ensuring that foreground contamination remains 
confined to outside the EoR window, safeguarding the 
potential of current generation experiments to make an 
exciting first detection of the EoR within the next few 
years. 
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